#################################################
####       Replication Code for              ####
####  A Re-Assessment of Reporting Bias in   ####
#### Event-Based Violence Data with Respect  ####
####         to Cell Phone Coverage          ####
#################################################
####                                         ####  
####                1/18/2017                ####         
####    (2) Generate Parameter Estimates     ####  
####     Cross-Sectional Window Analysis     ####  
#################################################


library(haven)
library(foreign)
library(ggplot2)
library(plyr)



### for figure 1a)
### high estimate, weidmann replication
load("africa_shape_grid_windows_high.RData")
load("W_windows.Rdata")


variable_dv = names(africa_shape_grid@data)[231:275]
africa_shape_grid@data$spat_lag = as.numeric(mat%*% africa_shape_grid@data$conf2007_dum)


iterator = seq(1:45)
model_func= function(i){
  africa_shape_grid@data$var = africa_shape_grid@data[, variable_dv[i]]
  m1 = glm(var ~ pre2000_count+bdist1+capdist+pop2005+mnt+irri+gcppc00+cell_dum_07 + spat_lag, family = binomial(link  = "logit"), data = africa_shape_grid@data)
  coef = summary(m1)$coef[["cell_dum_07","Estimate"]]
  low = summary(m1)$coef[["cell_dum_07","Estimate"]] - qnorm(0.975)*summary(m1)$coef[["cell_dum_07","Std. Error"]]
  high = summary(m1)$coef[["cell_dum_07","Estimate"]] + qnorm(0.975)*summary(m1)$coef[["cell_dum_07","Std. Error"]]
  res = c(coef, low, high)
  return(res) 
}

results_Windows_high = lapply(iterator,FUN = model_func)
save(results_Windows_high,file="results_Windows_high.rda")



# Generate results for Models in results_SpatWindows.rda
##set your working directory that contains all the files

load("africa_shape_grid_windows_best.RData")
load("EventsPerWindow.rda")
load("W_windows.Rdata")


variable_dv = names(africa_shape_grid@data)[231:275]
spat_lag = names(africa_shape_grid@data)[141:185]

iterator = seq(1:45)
model_func= function(i){
  africa_shape_grid@data$var = africa_shape_grid@data[, variable_dv[i]]
  africa_shape_grid@data$spat_lag = as.numeric(mat%*% africa_shape_grid@data[,spat_lag[i]])
  m1 = glm(var ~ pre2000_count+bdist1+capdist+pop2005+mnt+irri+gcppc00+cell_dum_07 + spat_lag, family = binomial(link  = "logit"), data = africa_shape_grid@data)
  coef = summary(m1)$coef[["cell_dum_07","Estimate"]]
  low = summary(m1)$coef[["cell_dum_07","Estimate"]] - qnorm(0.975)*summary(m1)$coef[["cell_dum_07","Std. Error"]]
  high = summary(m1)$coef[["cell_dum_07","Estimate"]] + qnorm(0.975)*summary(m1)$coef[["cell_dum_07","Std. Error"]]
  res = c(coef, low, high)
  return(res) 
}

results_SpatWindows = lapply(iterator,FUN = model_func)
save(results_SpatWindows,file="results_SpatWindows.rda")


# generate results for 
load("new_win_high.rda")
load("EventsPerWindow2_high.rda")

country_means <- ddply(new,~cow,summarise,cmeans=mean(cell_dum_07))
new = join(new, country_means,by="cow")
### running the windows mdels

modelfunction = function(y){
  new$var = new[,y]
  #standard logit, country mean
  m1 = glm(var ~ bdist1+capdist+pop2005+mnt+irri+gcppc00+cell_dum_07 + cmeans, family = binomial(link  = "logit"), data = new)
  coef = summary(m1)$coef[["cell_dum_07","Estimate"]]
  low = summary(m1)$coef[["cell_dum_07","Estimate"]] - qnorm(0.975)*summary(m1)$coef[["cell_dum_07","Std. Error"]]
  high = summary(m1)$coef[["cell_dum_07","Estimate"]] + qnorm(0.975)*summary(m1)$coef[["cell_dum_07","Std. Error"]]
  res = c(coef, low, high)
  
  return(res) 
}

variables_dummy = names(new)[c(231:275)]

results_windows = lapply(variables_dummy,FUN = modelfunction)
save(results_windows,file="results_windows_cmean_high.rda")


##### windowed analysis with first order by best then high estimate to break ties
##load data
llCRS<-CRS("+proj=longlat +ellps=WGS84")
load("africa_shape_grid_windows.RData")
load("W_windows.Rdata")
conflict = read.csv("UCDP GED v.1.0-2011.csv")  

conf2008<-subset(conflict, Year==2008)

windowsStart <- round(seq(1,round(dim(conf2008)[1]/2),length.out = 45))
windowsEnd <- round(seq(round(dim(conf2008)[1]/2),dim(conf2008)[1],length.out = 45))
africa_shape_grid@data$spat_lag = as.numeric(mat%*% africa_shape_grid@data$conf2007_dum)

conf2008 <- conf2008[order(conf2008$Best_est, conf2008$High_est),]
iterator = seq(1:45)

model_func= function(i){
  sub = conf2008[c(windowsStart[i]:windowsEnd[i]),]
  coordinates(sub)=~Lon+Lat
  proj4string(sub)=llCRS
  EventsInWindow = mean(sub$Best_est)
  sub <-SpatialPoints(sub,proj4string=llCRS)
  points_sub2008 <-over(africa_shape_grid,sub,returnList=TRUE)
  for (j in 1:10674){
    africa_shape_grid@data$conf_window[j] <- ifelse(length(points_sub2008[[j]])>0,1,0)
  }
  mod1 = glm(conf_window ~ pre2000_count+bdist1+capdist+pop2005+mnt+irri+gcppc00+cell_dum_07 + spat_lag,data=africa_shape_grid@data, family = binomial(link= "logit"))
  coef = summary(mod1)$coef[["cell_dum_07","Estimate"]]
  low = summary(mod1)$coef[["cell_dum_07","Estimate"]] - qnorm(0.975)*summary(mod1)$coef[["cell_dum_07","Std. Error"]]
  high = summary(mod1)$coef[["cell_dum_07","Estimate"]] + qnorm(0.975)*summary(mod1)$coef[["cell_dum_07","Std. Error"]]
  res = c(coef, low, high)
  return(res)
}
results_SpatWindows_BestHigh = lapply(iterator,FUN = model_func)
save(results_SpatWindows_BestHigh,file="results_SpatWindows_BestHigh.rda")


